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We introduce a method where successive coordinate transformations are applied to decrease the 
error in the adiabatic master equation resulting from truncation in the local adiabatic parameter. 
Our method reduces the non-physical behaviour stemming from the lack of complete positivity. The 
strong environment-induced relaxation at high Cooper pair pumping frequencies leads to adiabatic 
ground-state pumping only in the lowest-order approximation. We illustrate the robustness of the 
frequency where the adiabaticity breaks down using the high-order theory and show the emergence 
of an optimal environmental coupling strength, for which ideal pumping is preserved for the highest 
frequency. Finally, we study the effect of quantum interference on the pumped current and give an 
estimate for the relaxation rate of an experimentally measured system. 
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I. INTRODUCTION 

Detection and manipulation of geometric phases^"— in 
superconducting quantum devices has been an area of ac- 
tive research in recent years with one of the ultimate goals 
being the ability to realize holonomic quantum gates^. 
Even though alternative methods to experimentally gen- 
erate and detect the geometric phases in such devices 
have been proposed and realizedi^, the link they have to 
Cooper pair pumping has attracted major interesii^rJ^ as 
it provides means of detecting the phases by a measure- 
ment of the dynamically and geometrically transferred 
charges. A measurement scheme based on one supercon- 
ducting island with two tunable Josephson junctions, the 
Cooper pair sluice in a superconducting loop, has been 
introducedi^ and experimentally realized^S. Recently, a 
system based on a similar device structure has been pro- 
posed to execute fully geometric quantum computing us- 
ing non-Abelian phases^ii^. 

Even though the principles of operation for Cooper 
pair pumps are well-known, accounting for system- 
environment interactions has been a work in progress. 
Cooper pair pumping being essentially a coherent pro- 
cess, a proper description of the operation must include 
the effect of the external environment. The usual meth- 
ods for describing the dynamics of open steered quantum 
systems^^""— have been shown to violate the conservation 
of physical observables such as the electric charge and a 
variety of means have been employed in an attempt to 
properly describe the dynamics^^"— . Recently, it was dis- 
covered that for a consistent description using the master 
equation for the reduced density matrix of the system, all 
the non-secular terms must be included to enable relax- 
ation to a proper basis and to ensure conservation of the 
pumped charge in Cooper pair pumping^li^. The same 
master equation has also been derived using superadia- 
batic basest. 

In this paper, we employ the recent methodology of 
superadiabatic bases^^ in deriving the master equation 
for steered systems and apply it after multiple coordi- 



nate transformations of the time-local basis. This results 
in a master equation where the truncation error related 
to the local adiabatic parameter is potentially decreased 
as a function of the number of coordinate transforma- 
tions. We apply our method to the problem of Cooper 
pair pumping. We show that for the zero-temperature 
environment and fast pumping, the ground-state adia- 
batic evolution revives in the relaxation dominated re- 
gion only using the adiabatic basis, that is, in the lowest 
order of our description. Furthermore, we show that the 
overestimation of the pumped charge caused by the non- 
positivity of the density matrix, is alleviated dramati- 
cally by the usage of high-order bases. We simulate the 
breakdown of adiabaticity with increasing pumping fre- 
quency and show the emergence of an optimal coupling 
strength preserving ideal pumping up to the highest fre- 
quency. We present a condition for the highest transition 
probability caused by constructive interference between 
driving-induced excitations generated at different times 
and show that it corresponds to the downward resonance 
peaks in the pumped current. Finally, we obtain an es- 
timate for the relaxation rate of the device employed in 
the experiments of Ref. [13 to pump Cooper pairs. 

The structure of this paper is as follows. In the next 
section, we introduce the model describing a driven quan- 
tum system and demonstrate our method of defining suc- 
cessive effective Hamiltonians by coordinate transforma- 
tions. In Sec. mil we write a master equation for the ma- 
trix elements of the reduced density matrix of the system 
taken in an n-times transformed time-dependent basis. 
In Sec. IIVI we use the master equation to model Cooper 
pair pumping. Furthermore, we simulate previous exper- 
iments on the pumped current and derive an estimate for 
the relaxation rate of a measured superconducting sys- 
tem. We conclude the paper in Sec. [Vj 



II. MODEL 

We study a quantum system with a Hamiltonian Hs 
which depends on a set of real control parameters {qk} 
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that vary in time. The system is assumed to be interact- 
ing with its environment such that the total Hamihonian 
is 

H{t)^Hs{t) + Vit)+HE, (1) 

where V{t) is the coupUng between the system and its 
environment and He is the Hamihonian of the environ- 
ment. We assume that the couphng is of the generic form 
V = A^X{t), where A is the system part of the couphng 
operator and X{t) acts in the Hilbert space of the envi- 
ronment. Let |m; q{t)) be the instantaneous eigenstate of 
Hs{t) and Em{t) the corresponding eigenenergy defined 
by Hs[q{t)\ \'m\q{t)) = E^[q{t)] \m]q{t)). In the context 
of adiabatic evolution, {\m;q{t))} is referred to as the 
adiabatic basis. We assume that the adiabatic states are 
normalized and non-degenerate. 

Let the Hamihonian Hs{t) be diagonalized in a fixed 
time- independent basis using the eigendecom- 

position as Hg'\t) = D\{t)Hs{t)Di{t), implying that 

{n f\H^g\t)\m f) ^ E„i{t)Snm- We define a similar trans- 
formation for the total density operator p{t) in the 
Schrodinger picture as p^'^\t) = b\{t)p{t)bi{i). It fol- 
lows from the Schrodinger equation that the evolution of 
p^^'^it) is governed by the effective Hamiltonian for the 
adiabatic basis 

H'^^'>{t)^H^s'\t) + f'Mt)+V^^\t)+HE, (2) 

where V'^^^t) = t)\{t)y {t)t)i{t) = t)\{t)ADi{t) ® X{t) 

and wi{t) = -ib\{t)t)i{t). The eigenbasis of H^s\t) + 
hwi (t) is usually referred to as the superadiabatic basis. 

We can further define a unitary transformation D2{t) 
making Hg'\t) + hwi{t) diagonal in the fixed basis^^. 
Thus the evolution of the density matrix p'^' = 
b\p^^^b2 = b\b\pbib2 is governed by the effective 
Hamiltonian for the first superadiabatic basis 

^(2) {t) = kf (t) + hw2it) + V^^^ (t) +He, (3) 

where H^g\t) = Dl{t)[H^g\t) + hwi{t)]b2{t), h^\t) = 

I)l{t)V^^\t)b2it), and W2 = -ib\{t)b2{t). This 
method of successive coordinate transformations can be 
continued to yield for the (n — l)th superadiabatic basis 
an effective Hamiltonian of 

= + hwn + + He, (4) 

where H^ = bl[H^s~^'^ + hwn-i]b^, V^^") 

(11^2 A:)^^IV^i^i(nr=2 A) and Wn = -iDlA, where 
we omitted explicitly marking the temporal dependence 
of the operators for clarity. The operator product is 
defined as n"=2 ~ £^2-03 • • ■ -D„_i-D„. If we de- 
fine i)^"^ = 117=2 D^ for n > 2 and i)^"^ / for 



n = 1, the density operators governed by the Hamiltoni- 
ans in Eqs. ([2]), (jS]) and (jlj obtain a more universal form 

— [b'^g^)'^ b\pbib^g \ Defining successive diagonal- 
izations in this manner proves useful in Sec. IIIII as the 
recently derived master equation^"— can be applied to 
solve the system dynamics using these high-order bases. 

The iterative method described here is an adaptation of 
Berry's concept^! he later referred to as adiabatic renor- 
malization^. It is based on the idea that each transfor- 
mation rotates the basis we use to describe the system 
dynamics ever closer to the exact evolving closed system 
state, that is. the time-dependence of the transformed 
system Hamiltonian is suppressed after each rotation. 
After n transformations, we define the time-dependent 
basis as {bib'^g^ I™/)}- This approach generally works 
only in the restricted sense, that is, after a number of 
iterations, the following rotations will not allow us to de- 
scribe the dynamics of the system more accurately^. 

Finally, we introduce the local adiabatic parameter as 
ai{t) = h\\wi{t)\\/ /S.{t), where we compare the Hilbert- 
Schmidt norm of the operator arising from the adiabatic 
evolution \ \wi{t)\\ = ^JTYs{wl{tywl{t)} to an instanta- 
neous minimum energy gap in the spectrum A(t). Here 
Trg denotes the trace over the system degrees of freedom 
and in the following we will use Tr^; to denote the trace 
over the environment degrees of freedom. The parameter 
ai{t) should give a good estimate for the degree of adi- 
abaticity of the evolutio n'^^i'^^ . In cyclic evolution with 
the period T, the parameter scales as 1/T and, thus, in 
adiabatic evolution we should have ai{t) <^ 1. 



III. MASTER EQUATION 

We consider an adiabatically steered two-level quan- 
tum system weakly coupled to its environment. We 
denote the ground and excited states of Hs in the 
Schrodinger picture as \g) and |e), respectively, with cor- 
responding eigenenergies Eg and E^- Using the interac- 
tion picture approach, a master equation was derived to 
describe the dynamics of such a system in Refs. [s^^lll up 
to the linear order in ai{t) and the quadratic order in the 
system-environment coupling. The method of derivation 
employed in Ref. |3^ is our starting point for developing 
a numerical scheme for obtaining a more accurate de- 
scription of the dissipative system dynamics. In Ref. l36l . 
a master equation for nonsteered systems was used in 
conjunction with the effective Hamiltonian in Eq. ([S]) to 
derive the leading-order master equation under steering. 
However, a similar derivation can be carried out using 

^(") for any n to obtain a master equation for the matrix 
elements of p'"' . Notice that even though the method of 
defining successive coordinate transformations can be ap- 
plied to a system with arbitrary number of energy levels, 
we constrain ourselves to the two-level case. This is prac- 
tical since our main goal is to explore the implications of 
applying our scheme compared to previous results^ii^^. 
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We define the reduced density operator of the system 
as pg"^ = Tr£;{/5("^} so that its diagonal element becomes 
Pas' — {^\P^s\^) ^"^^ the off-diagonal element p'g]) ~ 
(0|pg'''|l), where {|m/ = 0) , |m/ = 1)} is the relevant 
fixed basis. These are simply the matrix elements of the 
usual density operator of the system in the Schrodinger 
picture taken in a time-dependent basis {15^"^) , |e'"^)}, 
where = DiD^g^ |0) and |e(")) = DiD^g^ This 

is the rotated basis obtained through the iterative pro- 
cedure we described in Sec. |IT1 We emphasize that the 
basis states are not obtained using a perturbative ex- 
pansion in the local adiabatic parameter and, thus, each 
iteration generally alters them by terms of all orders of 
al{i^r^^^. However, if the method described in Ref. [36l 



is applied in the n-times transformed basis, the error in 
the resulting master equation is defined by a perturbative 
expansion in Wn- The norm of Wn decreases, in the re- 
stricted sense, with increasing n as the time-dependence 
of the transformed system Hamiltonian is suppressed. In 
regard to depicting the actual dynamics of the system, 
the only issue relevant to the selection of the basis is that 
the time-evolution of the density matrix elements can be 
accurately described using it. Thus, we can exploit the 
nth iterative basis and define a master equation up to the 
quadratic order in the system-environment coupling and 
to the first order in Q;„(i) = ft||w„(i)||/cjQ"''(<), where 

^(n) _ ^(n) ^ ^(n) ^^^^ ^^^^ ^(n) ^ and 

4"' = (0|^i"^|0>, as 



+ 5ie(m("))Ke(p(:))]5(0)m("' - 2 ^^^°^ ""^^ ^^""^ ^ {[5m(m^"))Sm(z.ff) + Ke(m("')5Re(u;(^))] 

"01 

X [3m(4"')3m(p(;^)) + 5Re(m("))Ke(p(^))]} + ^ ^^W - ) - S{J-^) |c^^(^(»))c^^(^(n)) 



"i 



(n) 
01 



"01 



and 
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2S'(0) - S{-uj, 



ff)-^(4l\ ^W)2 r„) r ) ^ ^ ^(0)-g(4l\ ^(n)).^(n) 



(n) ^1 ' Se '^SS ' (n) 

"oi "oi 



_ "oi) M"oi) |c^^^(„4n))j^^(^(»)) _ 3„,(^(^))5Re(mr))} 

"01 



01 

I 
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3m(m("')3m(ii;(;')) + 5fte(m^"')5?e(i«(:))]p(^)}. 



(6) 



Furthermore, we denote rri\^ = (g'-"-' | Ajg'-"'), 

m^"' = (g(")|i|e(")), = -z (0|I)tA|0), = 

and w^gi) = -i (0|i)J,Ai|l)- The re- 
duced spectral density of the noise source is defined as 
S{uj) = dTTrE{pEX{T)X{0)}e''^yh^. Similarly to 
Refs. (3343 6^ . we assume that the system is in the Markov 



regime, the system time scales are longer than the envi- 
ronment autocorrelation time leading to neglecting the 
Lamb shift, and the approximation of adiabatic rates ap- 
plies. These assumptions and the time scale separation 
they lead to are described in detail in Ref. [s^ 

As described above, the benefit of defining the suc- 
cessive coordinate transformations of the Hamiltonian is 



4 



that the corresponding master equation is up to the first 
order in Q;„(t), thus describing the evolution of the system 
more accurately, in the restricted sense, as n increases. 
Defining a master equation of arbitrary order in ai (t) 
using the original methods^"— is possible, but the effort 
required renders such derivations highly unpractical. 

It has been shown^ii^ that assuming a zero- 
temperature environment and taking the quasi- 
stationary limit, the master equation in the lowest order 
-2^andp(i'- ~-^''> '-^'^ ' r,<.2. 



leads to p"gg 



1- 



.(2) 



•Jge I W| 

n(2) 



01 



This translates to p'gg = 1 + O(a^) and pfe =0-1- 0(ai) 
showing that in the first order in ai, the density 
matrix fP'^ describes the evolution of a pure state. 
This is a remarkable result validating that the master 
equation in Refs. (33 - [36| ensures relaxation to Lg*-^-*) 
up to the first order in a\. Similarly, using our master 
equation for p*^") ensures that the relaxation takes 
the system to l.g'"'*'^-') up to the first order in a„. 
Especially, in the limit n — ^ c», the rotational terms 



,(") 



, k,l € {g, e}, in the master equation vanish and the 



basis , |e("))}|. 



fully describes the steering 



assuming that the process of basis rotations converges. 
The requirements for the convergence or the number 
of transformations up to which the iterative procedure 
suppresses the time-dependence of the Hamiltonian of 
the systen*i are not studied in this paper, as it turns 
out in Sec. IIVI that a small number of transformations 
allows one to capture the main effect of this scheme. 



IV. COOPER PAIR SLUICE 
A. Definitions 

We introduce the Cooper pair sluicei^ as a physical 
realization of a steered two-level system. The charge 
pumped through the sluice establishes a connection to ge- 
ometric phases^iiiiiiiii^i^S acquired during the adiabatic 
evolution and provides a physical observable. We aim 
to study and improve on the recent theoretical pumping 
result a^^i^^ using the high-order effective theory. 

The Cooper pair sluice shown in Fig. [IJa) is comprised 
of a superconducting island separated by two SQUIDs^, 
each involving two Josephson junctions. If we assume 
that the self-inductances of the SQUID loops are negligi- 
ble, the SQUIDs operate as tunable Josephson junctions 
whose Josephson energies Ejk{^k), where k G {l,r}, are 
determined by the external fluxes threading the loops. 
The sluice Hamiltonian is 



— -Eji cos U> - — 
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Hs = Ec{n~ngY -Ejr cos 

'(7) 

where the Coulomb energy for one excess Cooper pair 
is Ec = 2e^/CE and the gate charge is defined in units 
of 2e as Ug = VgCg/{2e). Here Cg is the gate capac- 
itance and Cs stands for the total capacitance of the 
island. The operator describing the phase on the island 



is (j) = {ifr — 'Pi)/'2 and its canonical conjugate, the op- 
erator describing the number of excess Cooper pairs on 
the island, is n = —id^. The gauge-invariant phase dif- 
ference over the device if = tpr + ifi is determined by 
ip = 2tt^ I $0 and kept constant during the evolution. 

We denote = —id^^ {k e {I, r}) as the Cooper pair 
number operator of the fcth SQUID and write the average 
value of the current through the fcth SQUID as^ 

{h) = ^ {Trs{ps[nk. Hs]} + Tr{p[nfc, V]}) , (8) 

where Ik denotes the respective current operator. If the 
environment does not directly induce a current, that is, 
[fik, t^] = 0, Eq. ([5]) defines the usual current operator— 



h. 



2ei 
IT 



2e dHs 
h dipk 



(9) 



If we set Ec ^ max{Eji, Ejr} and Ug w 1/2, the dy- 
namics are accurately described by the two lowest charge 
states allowing us to apply the preceding two-state the- 
ory. We denote |0) and |1) as the states with no and 
with one excess Cooper pair on the island defining our 
fixed basis. We study the capacitive coupling of the en- 
vironment to the system by introducing voltage fluctu- 
ations SVg{t) at the gate of the sluice^ii^. The cou- 
pling operator then becomes Vg = —eg&z €D ^Vgit) where 
az = |0) (0| — |1) (1| and g = Cg/Cj: denotes the strength 
of the coupling. The coupling operator has been selected 
traceless in the two-state basis by adding an operator 
comparable to the identity operator to adopt a conven- 
tion used in the derivation of the master equation^. Such 
a selection can be applied to any coupling operator and 
it does not reduce the generality of the master equa- 
tion. Since [7ifc,F] = 0, Eqs. ([7]) and ^ imply that 
ik = 2eEjkSin{ipk)/h. 

Assume that the noise source is a resistor in thermal 
equilibrium; a situation which can be engineered in the 
physical realization of the sluice. If we consider the volt- 
age noise of a resistor grounded at one end and connected 
to the gate by a low impedance circuit at the other end, 
the reduced spectral density of the noise source at the 
gate becomesii S'(cj) = 2i?ga;/[ft(l-e-'*"/(''«'^«))], where 
Rg is the effective resistance of the noise source and Tr 
is the resistor temperature. We note that the detailed 
balance condition S{uj) = e'^'^/(''«^«)S'(-a;) applies. Fur- 
thermore, we introduce dephasing to the system by as- 
suming that 5(0) = 2kBToRg/h^, where Tq is the effec- 
tive dephasing temperature. 

We denote the matrix elements of the current operator 
of the fcth SQUID by ^"^^ = (r(") |/fe|s(")). The expec- 
tation value of the current using the adiabatic basis is 

ih) = + P^I^L + 2Jfte(pW4:4), (10) 

since Ik is Hermitian. The flrst two terms are the dy- 
namic supercurrent through the junction and the third 
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FIG. 1. (Color online) (a) Circuit diagram of the Cooper pair sluice. The fluxes threading the left and right SQUIDs are 
denoted by $; and $r, and $ denotes the total flux threading the large superconducting loop. The phase differences over 
the SQUIDs are marked by ipi and i^r and they are defined in the direction specified by the arrows in the figure. The gate 
capacitance Cg is used to manipulate the island charge with the gate voltage Vg. (b) The time dependence of the control 
parameters Eji, Ejr, and Ug during a pumping cycle. The residual values ej, tr, and 77 allow for non-ideal SQUIDs and gate 
control. 



term describes the geometric part of the currenti2i2£ii2ii^. 
The pumped charge corresponding to the geometric con- 
tribution becomes 



iG,k 



23?e 



t+T 



(11) 



where T is the length of the closed cycle. 



The definition for the different terms in the average 
current only applies directly for the adiabatic basis. This 
implies that if we pursue to obtain the geometric current 
using higher order bases, should be written using the 
density matrix elements in the basis where the evolution 
takes place in our calculations. The adiabatic density 
matrix element can be written as 



fl2) 



Using this, we can rewrite the integrand in the pumped 
charge 



\D 



(") 



l)(l|(4"))t|l))} 



.5Re{p(:)}5Rc{/W^((0|i)^")|0)(l|(4"V|l) 
3m{p(^)}3m{/« f(0|I)("^|0)(I|(^("V|I) 



(13) 



-(0\D 



(n) 
S 

r(i) 



1) (01(4"')^))} 



+ KeK:4 (0|4"^|l)(l|(I?r)'|l)}- 

Equation enables us to calculate the pumped charge 
when the time-evolution of the nth basis density matrix 



(")At 



is known. Note that for the adiabatic basis, Eq. 
reduces to the form corresponding to Eq. pip . 



B. Effect of the environment on the pumped charge 

We use the parameter cycle shown in Fig. [TJb) for the 
pumping and ensure the smoothness of the parameter 
functions in time using trigonometric interpolation di- 
viding the total cycle time into 201 equidistant points. 
This is a necessary step for the simulations as exploiting 
the high-order bases requires that the high-order tempo- 
ral derivatives of the parameter functions stemming from 
Wn are non-divergent. The dynamics of the quantum sys- 
tem are solved numerically from Eqs. ([5]) and ^ utilizing 
the effective Hamiltonians introduced in Sec.|lTl The den- 
sity matrix and any physical observables are recorded in 
the steady state, that is, after sufficiently many cycles 
such that the system evolution of consecutive cycles is 
identical. 

We begin by studying the effect of using the higher- 
order bases when describing the dynamics of the sluice 
in a zero-temperature environment. Using the lowest- 
order approximation n = 1, the analytical result of the 
environment inducing ground-state evolution in the adia- 
batic limit has been demonstrated using numerical simu- 
lations of the pumped charge^^i^^. Additionally, increas- 
ing the pumping frequency has been shown to induce re- 
gions where either the nonadiabatic transitions or relax- 
ation dominates depending on the ratio oti/g. However, 
it turns out that a region in the (ai,5)-space emerges 
where the pumped charge is unphysically overestimated. 
This is a direct result of the master equation not strictly 
ensuring the positivity of the density matrix in any fi- 
nite order, that is, it does not reduce to the standard 
Lindblad form. 

We present the pumped charge using the basis 
{|g(")) Je("))} with n = 1,2,3 as a function of the 
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coupling strength g for different pumping frequencies in 
Fig- assuming a zero-temperature environment. We 
explore the regime where the strength of the environmen- 
tal coupling is small to ensure that we remain close to 
the weak coupling limit. In the adiabatic region (/ = 10 
MHz), all orders of approximation indicate ground-state 
pumping for all environmental coupling strengths as pre- 
dicted in Sec. Hill By increasing the pumping frequency, 
we observe the emergence of the two pumping regimes 
mentioned above. Additionally, Figure [H^a) illustrates 
how increasing the coupling strength does not lead to 
ground-state pumping beyond the adiabatic region for 
n > 1. The reason for this phenomenon stems from the 
structure of the superadiabatic bases. For nonadiabatic 
evolution, increasing the coupling strength leads to re- 
laxation to 

|^(n+l)^ up to the first order in a„ which 
translates to ideal pumping only for n = 1 as the sys- 
tem is forced to the solution of the adiabatic limit, from 
which the asymptotic solutions for n > I generally de- 
viate in all orders of the local adiabatic parameter. The 
adiabatic solution for the pumped charge assuming simi- 
lar SQUIDS e ^ Er = ei has been derived previouslyi^ as 
Qc = 2e(l — 2ecosiy9). To make the discrepancy between 
the high-order bases and the case n = 1 more visible, we 
present the pumped charge up to high coupling strengths 
in Fig.[2{b). The regime of increased coupling strength is 
beyond the range of validity of our approach but shows 
how the master equation properly displays the distinc- 
tions between the different bases as they approach the 
asymptotic relaxation dominated solutions. 

From Fig. [21 it is clear that the charge overestimation 
observed for n = 1 is alleviated by utilizing the high-order 
bases. As the high-order bases follow the exact evolv- 
ing closed system state more closely, the non-adiabatic 
transitions disturb the mixed state less. The major con- 
tribution in alleviating the lack of complete positivity is 
already given by n = 2 and the subsequent third basis 
rotation has little effect on the pumped charge in compar- 
ison. The effect of using the higher-order bases is evident 
from studying not only a set of observables but from the 
density matrix itself. We present the lower eigenvalue Kg 

of the steady state density operator p^"-* in Fig. [31 The 
non-positivity is greatly reduced and for some parameter 
values, completely removed by doing one or more further 
rotations beyond the adiabatic basis. 



C. Pumped current and the breakdown of 
adiabaticity 

So far, experimental results for Cooper pair pumping 
using the sluice have been scarcoi^iSSi^. However, the 
breakdown of adiabaticity has been observed by study- 
ing the pumped current Iq as a function of the pumping 
amplitude Grange = max{ng} — minjrig}, that is, at high 
pumping amplitudes the pumped current has been no- 
ticed to deviate from the analytical result in the adiabatic 
limit |/g| = 2encp/, where Ucp is the number of Cooper 



pairs transported ideally through the sluice per cycle^S. 
The number of transported Cooper pairs can be exper- 
imentally dictated by adjusting the gate voltage, more 
spesifically, by altering Grange so that the ideally trans- 
ported average current corresponds to \Iq\ = 2e7ii-angc/- 
Unfortunately, our master equation cannot be directly 
used to simulate the effect of altering the pumping am- 
plitude since it is defined in a two-state basis which re- 
quires that Ug remains approximately half-integer during 
the evolution. Any selection of the two charge states as 
the fixed basis enables the maximum pumping amplitude 
of one Cooper pair per cycle. 

Even though we cannot modify the pumping ampli- 
tude, we can still simulate the breakdown of the adia- 
baticity by altering the pumping frequency. If we assume 
that the deviation from the adiabatic behaviour in the ex- 
periments is due to increase in the pumping speed caused 
by the amplitude growth using a constant cycle time, the 
effect of decreasing the total cycle time should be similar. 
We choose a frequency range beyond the strict adiabatic 
limit to model experiments carried out with finite cycle 
times. The pumped current is shown in Fig. [H^a) using 
n = 1 and in Fig.lDJb) using n = 3. The behaviour of the 
pumped current with high frequencies is only suggestive 
but clearly indicates that the effect of the number of co- 
ordinate transformations increases with the pumping fre- 
quency. The physically most relevant features are found 
near the point where the adiabaticity breaks down. As 
we utilize a more accurate description of the dynamics, 
the point where the adiabaticity of the system is bro- 
ken becomes more robust against changes in the envi- 
ronment. Furthermore, we observe the emergence of an 
optimal coupling strength, with which the ideal ground- 
state pumping is conserved up to the highest frequency. 
This feature could have also been anticipated from Fig. [5] 
where we observe the emergence of a maximum pumped 
charge as a function of the coupling srength for any given 
frequency far from the adiabatic limit if n > 1. This 
corresponds to the coupling strength, below which the 
non-adiabatic transitions reduce the pumped charge and 
above which the relaxation takes the system away from 
the solution in the adiabatic limit. The optimal coupling 
strength can be probed and exploited utilizin g th e envi- 
ronment engineering scheme presented in Ref |35l 

We turn our attention to modeling the breakdown 
characteristics of the actual experimentally pumped cur- 
rent in Ref. [20l . The simulation is motivated by the 
above discussion on the pumping speed and we es- 
tablish an equality between the experimental pumping 
speed T^rango ^ /oxpj whcrc /oxp is maintained at a con- 
stant value, and the pumping speed in our simulation 
nf^^^ X f , where n?!™ ^ is a constant to ensure that the 

1 angc •' ' 1 angc 

two-state approximation holds. A comparison between 
the experimental results and our simulations is presented 
in Fig. [5l 

The system parameters used in the simulations are 
estimated from the experiments and can contain up to 
20% error. To simulate the experiments, the aforemen- 



7 




— n=l, /=10 MHz 
^n=l, /=50 MHz 
/=100 MHz 
f=10 MHz 
/=50 MHz 
-^n=2, / = 100 MHz 
f=W MHz 
/=50 MHz 
— n=3, / = 100 MHz 



0.12 0.14 




— n 


=1, /= 


100 


MHz 


-^n 


=1, /= 


300 


MHz 


-^n 


=2, /= 


100 


MHz 


-T-n 


=2, /= 


300 


MHz 


-^n 


=3, /= 


100 


MHz 




=3, /= 


300 


MHz 



am 0.05 



0.15 
9 



0.3 



FIG. 2. (Color online) Charge pumped through the Cooper pair sluice during one cycle assuming that the temperature of 
the environment is Tr = (a) with the pumping frequency / = 10, 50, and 100 MHz and (b) / = 100 and 300 MHz. The 
simulations are carried out setting the number of coordinate transformations to n = 1, 2, and 3. The physical parameters used 



in the simulations are To 
7] = 0.25 and (p = n/2. 
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FIG. 3. (Color online) Lower eigenvalue Kg of the reduced density operator after n = 1, 2, and 3 coordinate transformations 
from bottom to top (a) for / = 10 MHz and g = 0.01 and (b) for / = 100 MHz and g = 0.03. Inset of (a) shows Kg for / = 10 
MHz and n = 2 and 3 from bottom to top. Other physical parameters are identical to those used in Fig. [2] 
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FIG. 4. (Color online) Current pumped through the sluice during one cycle as a function of the pumping frequency / with 
coupling strengths g = 0.025, 0.05, 0.1 and 0.3 using (a) n =1 and (b) n = 3. Other physical parameters are identical to those 
used in Fig. (2] The linear solid line is the ideal pumping curve for the two-state model |7g| = 2e/. 
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FIG. 5. (Color online) (a)-(c) Pumped current per cycle given by the experiment (circles) and by the simulation (dashed lines). 
The linear solid line is the ideal pumping curve assuming that the number of geometrically transported Cooper pairs is given by 
the pumping amplitude \Ig \ = Zeji^^g;, /cxp . The parameters from the experiments are Ec/kB = 2 K, E'f^'^/kB = E^l^^ /kB = 



0.7151 K, er = ti^ 0.05, Lp = 7r/2, and /oxp = 14 MHz. For the simulations, we use Tr = 200 mK, Rg = 300 kf7, 



= 0.8, 



and r] = 0.25 giving n^"ngt, = 0.6. The effect of dephasing is studied by using (a) To = K and (b) To = 0.1 K, where the 
coupling strength is g = 0.025, 0.05 from bottom to top. In (c), we present the pumped current for To = 1 K and g = 0.025 
(solid), and the averaged pumped current for To = 0.1 K and g — 0.025 (dashed). The corresponding relaxation rate of the 
system (solid), and the averaged relaxation rate (dashed) are given in (d) as a function of time. Experimental data courtesy of 
M. Mottonen, J. J. Vartiainen, and J. P. Pekola. 



tioned equality defines the needed frequency by / = 
"range ^ 23.33 MHz SO that tlie point where the adia- 
baticity breaks down according to the experimental data 
would require a simulation frequency of approximately 
1.1 GHz. This is significantly higher than the breakdown 
frequency in our previous simulations and tests the limits 
of validity of our approach. 

Assuming both vanishing dephasing in Fig. [SJa) and 
non- vanishing dephasing in Fig. [SJb), the simulated 
pumped current exhibits strong oscillatory behavior al- 
lowing for weak predictability of the exact breakdown 
characteristics. The observed behavior is due to quantum 
interference between driving-induced excitations gener- 
ated at different times. The instantaneous transition 
probability is not only dependent on the energy gap but 
also on the phases accumulated during the quantum evo- 
lution. Especially, if the phase difference between two 
successive excitations in time is a multiple of 27r, the 
transition probability reaches its maximum due to con- 
structive interference and the geometric current obtains a 
downward resonance peak. Studying the phase accumu- 
lation during the quantum evolution allows us to estimate 



the resonance peak positions and compare them with the 
simulated pumped current. 

The time-evolution of the energy gap of the system is 
symmetric with respect to the mid-point of each pumping 
cycle corresponding to times t^^-^^ — {2i -|- 1) x T/2, i S 
{0, 1, 2 . . . }. In addition, the energy gap decreases during 
the two gate manipulations and as both manipulations 
in Figs.[5{a) and [Sib) are symmetric with respect to the 
degeneracy point Ug = 1/2, the energy gap reaches its 
minima at points we denote as ti = + x T/4. Hence, 
it is enough to study the phase accumulation between 
different ti to depict the resonance behavior. 

In adiabatic evolution, the number of coordinate trans- 
formations we perform has a profound effect on the ob- 
served accumulated phase, that is, each transformation 
takes us to a new basis in which the phase accumulates 
differently. More exactly, the closed system state |4'(i)) 
fohows the Schrodinger equation ih\ii{t)) = Hs{t) \^{t)) 
at all times. After m transformations, the state 
|*("0(t)) = Di{t)D'^g"'\t) |*(t)) follows the transformed 

equation i?i !*('")(<)) = (#^"^(t) -|- ^^^(t)) |*(™)(t)). 
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Assuming adiabaticity in this basis, that is, the exact 

evolving state remains in the fcth eigenstate of H^"^ (t) + 
hwmit), the m-times transformed evolving state can be 
written, similarly to Ref. [4^ as 

\¥^--Ht)) = exp (^^4^^\t) i £ drEi"^\r)^ |fc(™)(i)) , 

(14) 

where we assumed that |^'('"H*in)) = Ifc^^H^in)), 

(^^")(t) + ;iu;„(t))|fc(")(t)) = ^^'"^(t)|fc(™)(t)), and 

a^™^ (t) describes the geometric phase contribution. 

Using the familiar notation, this implies E^^^\t) ~ 

£;[."+^'(t) and |fc(™)(t)) ^ l)™+i(t) |fc/). The accumu- 
lated quantum phases can be obtained from Eq. ([Ti|) fol- 
lowing the derivation in Ref. |43| . However, in our dis- 
sipative calculations, after n transformations, we take 

hw^ as the perturbation. This means that for the 

dissipative simulations, the relevant reference frame af- 
ter n transformations is given by the eigenbasis of Hs {t) 

for n = 1 and Hg^~^\t) + fvWn-iit) for n > 1, that is, 

Et''\t) = Et\t) and \U--^\t)) = b^{t) \kf). 

Concentrating on the two-state model, we present the 
accumulated phases after n transformations. We as- 
sume that the transitions taking place at ti are instan- 
taneous and the system evolves adiabatically between 
them. The dynamically accumulated phase obtained dur- 
ing adiabatic evolution between ti and t^+i for the eigen- 
states in the relevant reference frame after n transfor- 
mations is given by Q^^\{ti,tij^i) = — J^^'^^ dtE'j^^ /h, 
k S {g,e}, and it is equal for any successive two points 
due to symmetry. The geometrically accumulated phase 

®G 1(^17 ^j-i-i) is not the usual Berry phase since ti with 
odd and even i describe different points in the Hamil- 
tonian space and, thus, the path traversed between ti 
and ti^i is not closed. Assuming adiabatic evolution, we 
can write the noncylic geometric phase acquired by the 
fcth eigenstate between two successive transitions using 
Eq. IHl) as^ii^ 

e^"i(t„t,+i) = arg{(fc("-i)(t,)|fc("-^Hi.+i))} 

+ 1 j dt {¥''-^\t)\dt\U''-^\t)) 

(15) 

= arg{{kf\Dl{t,)Dr,{t,+i)\kf)} 
+ 1 j^^^\t{kf\bi{i)b,,{t)\kf). 

The geometric phase defined in Eq. ([T5l) is gauge- 
invariant and, thus, only depends on the traversed path. 
Using our parameter cycle, the accumulated geometric 
phase for any successive two points is the same. Note that 
as the transformations suppress the time-dependence of 
the effective Hamiltonians, the geometric phase is de- 
creased as n increases. The difference in the total accu- 



mulated phase acquired by the eigenstates is 

Ae^'^(t,;,i,+l) 

= arg{(0|-Dt (^,;)A^(^.+l)|0)} - arg{(l|I)t (^^)A.(^^+l)|l)} 

+ j'"^ dt [z«(:)(t) - w^^^{t)] + dt J.'i^t). 

(16) 

For a full cycle, we have A0y' (ii, ^^+2) w 2{ip— 2e sin ip) + 

Jf'^^ ujQ^i(t), where we assume similar SQUIDs e = = 
ei, as the accumulated geometric phases become the 
Berry phases in the adiabatic basis^ii^. 

In addition to the phases accumulated in adiabatic 
evolution, we must account for any phase shifts occur- 
ring at the transition. We apply the adiabatic-impulse 
model^ in the relevant reference frame and describe the 
non-adiabatic transitions taking place at ti as instanta- 
neous processes. The model is based on the Landau- 
Zener approximation^ assuming that in the vicinity of 

each ti, Hsit) for ri = 1 and Hg" ^\t) + hwn-i{t) for 
n > 1 can be linearized in a fixed basis as 

H[^^(t) = -A„/2a, T Vnt/2a,, (17) 

where ctx = |0) (1| + |1) (0|, the upper sign corresponds 
to odd i and the lower sign corresponds to even i. 
The representation of the Hamiltonian exploiting a time- 
independent basis is usually referred to as diabatic. The 
tunneling amplitude A„ and the slope of the linearized 
energy bias w„ are assumed real. Our Hamiltonians do 
not exactly linearize into this form, but the off-diagonal 
elements obtain phase factors related to the phase differ- 
ence over the sluice for n = 1 and to the complex phase of 
Wge for n > 1. These phase factors can be accounted for 
by a transformation to a representation where the Hamil- 
tonian linearizes to the Landau-Zener form. The Landau- 
Zener transition probability from the ground state to the 
excited state is given by Pjj^ = exp(— 27r(5„), where the 
adiabaticity parameter is (5„ = A^/(4z;„). 

We can generally linearize the relevant Hamiltonian 
after n transformations in the diabatic representation in 
the vicinity of each ti as 

H(::)(t) = -A„(i,)e"^"(*')/2|0)(1| 

-A„(t0e-'^^"(*'V2|1)(0| (18) 
T v„(ti)t/2az, 

where 71 (t,) = aig{Ejr{t^)e-'^/^ + Eji{t,)e'^/^} and 

ln{U) = arg{— for n > 1. Additionally, 
we have Ai(t,) = \Ejr{ti)e-''^/'^ + Eji{ti)e'v/^\ and 

A„(ti) = \2h'wl^~'^\ti)\ for n> 1. Since the gate charge 
is altered symmetrically with respect to the degeneracy 
point Ug — 1/2 in Fig. [SI the linearization yields vi(ti) = 
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-2Ecng{U) and v„{U) = -2[E^;'-^\u)) + h4j'g-^\ti)\ 
for n > 1. Then the state of the system \i^{t)) evolves 

accordmg to ih\'4}{t)) = Hi^{t) \ip{t)) near ti. Defining 
a transformation 

U'^^(t) = e*^"W/2 |0) (0| + e-''^"(*)/2 ^ (19) 

yields that the transformed state — 

(C/(^)(0)t|^(t)) follows ^h\m ^ Ht^it)W)) 
near ti and we can apply the adiabatic impulse model. 
The evolution operator describing a single Landau-Zener 
transition is^ 

-^\e'^iiU)){sl:^iU)\ (20) 

+ v/l - Pll^e'-s^ I5[z^(i0)(3[z^(i0l, 

where At is a short time-step, (^^"'' = — 7r/2, 1^9 = 
7r/2+(5„(ln (5„ — l) + arg{r(l — F is the gamma func- 
tion, and {|(?L2 (^)) ' I^LZ (0)} is the instantaneous eigen- 
basis of i?Lz' (^) • "^^^ evolution operator describing the 
Landau-Zener transition in the original representation is 
obtained with a simple back-transformation iV("^(tj — 

At,u + At) if^^\u)Nl^^\u - At,u + At)iui;^\t,)y . 

Since the system is in the charging regime, di ^ 1 
implying that pjj^ <C 1. Furthermore, P^z^ decreases 
with increasing n since (5„ ^ c^n-i- Thus, with suf- 
ficiently large n, we are in the slow-passage limit^ 
and, additionally, the temporal suppression of the ef- 
fective Hamiltonians in the adiabatic renormalization 
yields that |g[z^(t)) |0) and \e'l2{t)) w |1) so that 
iV*^"^ {ti — At, ti + At) becomes approximately diagonal in 
the fixed basis. This implies that the transition causes 
an impulsive phase difference between the instantaneous 
eigenstates of the original effective Hamiltonian as 

= ^Llgit^ - t^ + ^t) " ^HeiU - At, t, + At) 

« — TT. 

(21) 

Note that in this limit, the transformation uj^iti) be- 
comes negligible and, hence, the phase factor 7„(t,;) does 
not appear in the impulsive phase difference. 

The reference times relevant to the interference are 
close to the transitions at ti ~ At and t^+i — At, and the 
time-evolution of the system is described by the evolution 
operator iV(")(i, - At,ti + At)U'^"-^t, + At,t,+i - At), 
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FIG. 6. (Color online) Comparison between the accumu- 
lated phase difference and the observed resonance peaks in 
the pumped current. The solid lines depict A0^' (ti, ti+i) 
for 71 = 1 and n = ?> from bottom to top. The dashed vertical 
lines are the resonance peak positions determined from the 
simulated pumped current in Fig. [Sjb) for g = 0.025. 



where (tin, ifin) describes the adiabatic evolution in 
the original representation. Hence, the phase differ- 
ence between the excitations is the sum of the impul- 
sive and adiabatic phase differences accumulated between 
the states Dn{t) |0) and Dn{t) |1)- This implies that the 
condition for the maximum constructive interference be- 
tween the successive excitations resulting in the maxi- 
mum transition probability is achieved with a frequency 
corresponding to 

^©LZ - + 

+ Ae^"^ (t, + At, t,+i - At) (22) 
= 2'kN, 

where TV e Z. Note that Ae^"^(ti + At,t,+i - At) = 

AQip\t„t,+,). 

In the limit of sufficiently large n, Eqs. (PTjl and ([^^ 
yield an approximate condition for the frequencies corre- 
sponding to the downwards resonance peaks as 

Ae^"^(t„t,+i) «2^ . (23) 

Note that AO^'' (ti, t^+i) is dependent on both the driv- 
ing frequency and the used basis. The oscillations should 
only be present with high frequency and low environ- 
mental coupling strength as decoherence destroys any 
low-amplitude interference effects. We present the ac- 
cumulated phase difference AQ^^\ti,ti^i) and compare 
it with the observed peak positions in Fig. [51 The reso- 
nance frequencies determined from the pumped current 
show an excellent agreement with the condition given in 
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Eq. ()23p . This reaffirms our assumption on the origin 
of the oscihations as being caused by quantum interfer- 
ence between excitations and offers a way to predict the 
resonance behaviour in possible experiments. Using the 
adiabatic basis also allows for a good estimate for low 
frequencies. Furthermore, we conducted simulations sim- 
ilar to the ones in Fig. [5] using (p = and obtained new 
downward resonance peak positions. Recalculating the 
accumulated phase difference allowed us to verify that 
the approximate condition in Eq. (j23p still applies with 
excellent accuracy. Further results regarding quantum 
interference in the Cooper pair sluice, especially with re- 
spect to applications in phase interferometry, have been 
derived recently^. 



Figures m a) and[5ljb). show that including dephasing 
decreases the amplitude of the oscillations. From the 
point-of-view of interference effects, this was to be ex- 
pected as finite dephasing time implies the excitations to 
lose some of their phase coherence during the evolution. 
The current variation in the experimental data is signifi- 
cantly lower implying that the effective dephasing rate is 
most likely higher in the experimental setup than in our 
simulation. However, to account for the large variation in 
the gate voltage, a model accounting for a higher number 
of energy levels would need to be introduced resulting in 
more complex transitions. Thus, we refrain from making 
further comparisons and study the breaking point of the 
adiabaticity. 



In an effort to suppress the oscillations, we perform 
the simulation for three different intervals for the gate 
charge maintaining Grange = 0.6 and take an average over 
the results. Moving the midpoint of the gate voltage in 
our control cycle slightly from the degeneracy point al- 
ters the temporal dependence of the energy gap, and, 
thus changes the positions of the resonance peaks. Tak- 
ing an average should then decrease the oscillations. We 
select the ranges as Ug £ [0.2 -I- 6ng;0.8 + Sug], where 
Sug € {—0.1,0,0.1}, so that the largest possible range 
for the gate variation is obtained while maintaining the 
two-state approximation with a reasonable accuracy. The 
result of such averaging procedure for g = 0.025 with de- 
phasing is shown in Fig. ^c). In Fig. ^c), we also give 
the pumped current with the same coupling strength and 
increased effective dephasing temperature to suppress the 
interference effects. Comparison with the experimen- 
tal data shows that we obtain a decent estimate for the 
breakdown characteristics with this environmental cou- 
pling strength using both methods. The temporal depen- 
dence of the corresponding relaxation rate of the system 

Teg = ITO2"'' P'^C'^oi'') ^^'^ °f relaxation rate aver- 
aged over the three different simulations is presented in 
Fig.EI^d) at / = 14 MHz. It is evident that the relaxation 
is strongest during gate operations due to non-adiabatic 
transitions, as expected. 



V. CONCLUSIONS 



We introduced and demonstrated a method of apply- 
ing successive coordinate transformations to describe ac- 
curately the dissipative dynamics of a steered two-level 
quantum system. Our method utilizes superadiabatic 
bases and the theory for nonsteered systems to obtain a 
master equation where the error resulting from the trun- 
cation of the perturbative expansion in the local adia- 
batic parameter is decreased. 

We applied our method to Cooper pair pumping and 
showed that in the adiabatic limit, all orders of approx- 
imation return the ideal pumping result. In the zero- 
temperature limit and increased pumping frequency, in- 
creasing the strength of the environmental coupling was 
shown to induce ideal ground-state pumping only in the 
lowest order of our description. Furthermore, using high- 
order bases was shown to reduce the overestimation of the 
pumped charge stemming from the non-positivity of the 
reduced density matrix of the system. This is due to the 
high-order bases tracking the exact evolving state more 
closely. The major effect of these corrections was shown 
to be captured by the basis obtained with two transfor- 
mations, i.e., the first superadiabatic basis. 

We studied the breakdown of adiabaticity by simulat- 
ing the pumped current with increasing pumping fre- 
quency. The high-order theory was shown to provide 
a more accurate picture of the robustness of the break- 
down frequency against changes in the environment. An 
optimal strength of the environmental coupling was dis- 
covered preserving the adiabaticity of the system for the 
highest pumping frequency. The recently proposed^ en- 
vironment engineering scheme can potentially be used to 
probe and exploit this optimal point. 

Finally, we applied our theory to model experimental 
pumping results similar to those of Ref. We altered 
the pumping frequency to simulate the increased pump- 
ing amplitude and observed oscillatory behaviour of the 
pumped current caused by quantum interference between 
driving-induced excitations generated at different times. 
We presented a condition for the highest excitation prob- 
ability due to constructive interference and showed that 
the observed downward resonance peaks in the pumped 
current accurately corresponded to this condition. Us- 
ing an averaging procedure and increasing the effective 
dephasing rate were methods that enabled us to finally 
present an estimate for the relaxation rate of the sys- 
tem. However, a many-state theory should be developed 
to facilitate more accurate predictions allowing us to al- 
ter the pumping amplitude. A seemingly valid approach 
could exploit the cyclic nature of the steering utilizing 
the Floquet theory, which would possibly allow for gen- 
eralizations beyond the adiabatic evolutional^. 
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